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We present a specially designed evolutionary algorithm for the prediction of surface reconstructions. This new 
technique allows one to automatically explore all the low-energy configurations with variable surface atoms and 
variable surface unit cells through the whole chemical potential range. Using this method to study the recon- 
structions of the GaN (1011) surface with and without absorbate oxygen, we find a non-intuitive reconstruction 
based on N3-trimers. The method can be generally applied to binary and multi-component systems. 



INTRODUCTION 

At crystalline surfaces, especially in semiconductors, atoms 
rearrange to form reconstructed configurations, which are dif- 
ferent from those in the bulk crystals. Reconstructions play a 
key role in surface properties, such as epitaxial crystal growth, 
catalysis, etc. Determining the reconstruction of semiconduc- 
tor surfaces is a long-standing problem. Despite progress in 
experimental techniques, it is still hard to obtain the atomic- 
scale details. Diffraction and microscopy can yield the peri- 
odicity and symmetry (TIE). However, the sampling space is 
astronomically large, and solutions may defy intuition, even 
for monoatomic systems: consider Si(lll) 7x7 BHU. For bi- 
nary and more complex cases, a variable stoichiometry at the 
surface adds a huge extra dimension. 

From the computational perspective, the main difficulty lies 
in that the number of candidate structures is extremely large. 
It is not possible to enumerate all the reasonable reconstruc- 
tion models just from chemical intuition. There have been 
a few efforts towards automatic prediction of the most sta- 
ble configuration for either a given stoichiometry |6-8], or a 
variable composition at constant chemical potential | 9 |. With 
these structural models available, the thermodynamical phase 
diagram can thus been constructed ifTOl . However, unbiased 
predicting reconstructions with variable stoichiometry in the 
whole chemical potential range has not been attempted to our 
knowledge. In this letter, we propose a new method which al- 
lows automatic exploration of all low-energy configurations, 
with variable surface atoms and variable reconstruction cells 
through the whole chemical potential range. We apply this 
method to GaN (1011) surfaces with and without added oxy- 
gen. 

THEORETICAL BACKGROUND 
Surface Energy and Chemical Potential 

Relative stability among candidate surfaces is determined 
by the formation energy ^formation 

^formation = ^tot — ^ref ~ / . n if^i, (1) 



where E tot and E YQ f is the total energy of the surface under 
consideration and of the reference cleaved surface, rii and jii 
are the number and chemical potential for each species, re- 
spectively. The chemical potential is the energy to add or re- 
move one atom from the system, assuming there is a reser- 
voir for each species to equilibrate with. The chemical poten- 
tials must satisfy constraints. For example in a simple binary 
AB compound, the chemical potential of A (or B) is bounded 
above by the chemical potential of elemental A (or B), 

Ma < Ea n , 
Mb<£b- {) 

The A and B atoms are in equilibrium with the substrate, 

Ma + Mb = ^ab- (3) 

These conditions lead to the range of the A component chem- 
ical potential to be 

E A b - E B < /i A < E A . (4) 
Therefore, Eq. ([T]) can be rewritten as only /i A dependent 

^formation = ^tot — ^ref — ^B^AB ~ Ma(^A ~ ^b)- (5) 

The surface phase diagram shows which reconstruction is 
favored under a specific environment, i.e., temperature and 
chemical potential or partial pressure of each component. 
Each of these environmental conditions represents an inde- 
pendent dimension in the phase diagram. 

Evolutionary algorithm 

We employ an evolutionary algorithm (EA) to explore the 
surface structures. The algorithm is implemented in the US- 
PEX package, which has been successfully applied to various 
bulk materials lfTTHT4l . We use a representation in which the 
system has three parts: vacuum, surface and substrate. Vac- 
uum and substrate regions are pre-specified, while surface re- 
gion is optimized by the EA. The number of surface atoms 
varies up to the maximum for a given thickness of surface. To 
allow automatic exploration of reconstructions with different 
surface unit cells, the cell size is also variable. 
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FIG. 1 . Surface representation used in the evolutionary algorithm. 



Our algorithm is constructed in the following way. It first 
generates random structures with only the surface region be- 
ing varied. These initial structures are then relaxed, and 
ranked by its fitness, based on the energy. Structures with bet- 
ter fitness are more likely to be selected as parents to generate 
new child structures. We considered three ways to produce 
offspring. 

(i) Heredity: Two structures are chosen from the previous 
generation. They are randomly sliced at the same position in 
the surface unit cell. Then pieces from both parent structures 
are combined to generate the offspring. 

(ii) Mutation: One structure is chosen as parent. Similar 
to the implementation of softmutation in bulk crystals [15], 
the surface atoms are mutated according to the softest surface 
vibrational modes based on the bond-hardness model fT6l[T7ll . 
If the structure cannot be softmutated anymore, we switch to 
the coordinate mutation [ 15 ], which is to randomly mutate the 
surface atoms from a zero-mean Gaussian distribution. 

(iii) Transmutation: One structure is chosen as parent. 
Some atoms are transmuted to another type of chemical iden- 
tity. 

The offsprings, together with a few best structures from the 
previous generation, comprise the next generation. This pro- 
cess is repeated until no new ground-state structures are pro- 
duced for sufficiently many generations. As a test, we applied 
the EA search to study the known 2x1 reconstructions of dia- 
mond (100) and (111) surfaces. We tried 2 and 6 carbon atoms 
on a 2 x 1 surface cell. During the search, each structure was 
first relaxed with the Brenner potential ifTSl as implemented 
in the GULP code (T9J, and then more accurate energy evalua- 
tion was performed at ab initio level with PBE functional l20l 
as implemented in the VASP code l2T1l . Fig. [2] shows the best 
structures found in the search. The ground state of diamond 
(100) surface contains characteristic C dimers (Fig. [2^), while 
the most stable configuration of diamond (111) surface indi- 
cates a rather complex reconstruction. Different from the bulk 



structure, the surface atoms form the 5+7 ring topology from 
the side view (Fig. [2]}, it is noteworthy that such peculliar 5+7 
ring topology can even form three dimensional structure via 
cold compression of graphite (221 123)). Our results are in ex- 
cellent agreement with those shown in the previous literature 
l24ll . We also performed the predictions for 4 and 12 carbon 
atoms on a 2x2 cell, but only with the Brenner potential. Al- 
though Brenner potential gives wrong energy ranking, the true 
ground states were also found as a low-energy configurations 
in the search. With this encouragement, we proceed to the 
case of more than one specie and variable stoichiometry. 

Fitness function 

Our fitness function needs to indicate the relative stability 
of the structures with various surface stoichiometry and re- 
construction cell sizes. We illustrate this in Fig. [3] First, 
for two given surface configurations (I and II) which have 
different numbers of atoms on the same substrate cell, their 
relative energy difference depends on the chemical potential 
Ha according to Eq.( [5]). As shown in Fig. [3^, surfaces I 
and II should coexist at ha = Meq- Therefore, the stability 
regime for the two surfaces can be fully established. Surface 
I is stable when Hmm < Ha < Meq, while surface II is sta- 
ble when /i eq < ha < Mmax- For any other unstable con- 
figuration like surface III, the fitness can be viewed as the 
minimum energy difference comparing with the stable con- 
figuration, at all possible chemical potentials. The minimum 
condition is reached when ha = Heq- F° r practical imple- 
mentation it is useful to express this algebraically. We define 
E = E tot — E re f — n^E ab as the h a -independent term in 
Eq.([5]). This approach was introduced by Qian et al |25 ], and 
is widely used l26l . Versions (a) and (b) of Fig. [3] contain 
equivalent information. The stable structures appearing in the 
phase diagram form a convex hull. The slope of each section 
of the convex hull is either the boundary chemical potential or 
the equilibrium chemical potential where stable structures can 
coexist. We then choose the fitness of a surface structure to 
be its distance to the convex hull. The EA search then aims to 
optimize the convex hull, similar to our previously proposed 
approach of variable composition prediction for bulk crystals 
(121 E21- When different sizes of surface cell (m x n) occur, 
the scale factor AT cell = m x n should be taken into account 
and energies have to be normalized. 

RESULTS AND DISCUSSIONS 

We studied the reconstructions on the semipolar GaN 
(1011) surfaces in the presence of oxygen. As a global op- 
timization method, an EA search requires hundreds or thou- 
sands of individual structural relaxations. In our calculation, 
relaxations were done using DFT within the generalized gra- 
dient approximation (GGA) [ 20 ] using the all-electron projec- 
tor augmented wave (PAW) |28 ] method as implemented in 
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FIG. 2. Schematic image of (a) diamond (100) 2x1 reconstruction (red dashed line: face of conventional cube; white rectangle: 2x1 surface 
unit cell); (b) diamond (111)2x1 reconstruction (red rectangle: 2x1 surface unit cell). 
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FIG. 3. Illustration of the fitness function used in the evolutionary al- 
gorithm, a). Phase diagram as a function of /xa- b). Phase diagram as 
a function of (ua — tib). The vertices of the convex hull are the stable 
structures appearing in the phase diagram. The slope of each section 
of the convex hull is either the boundary chemical potential or the 
equilibrium chemical potential where stable structures can coexist. 



the VASP code [21 ]. We used the kinetic energy cutoff of 550 
eV for the plane-wave basis set and Brillouin zone sampling 
resolution of 2tt x 0.08 A -1 , which showed excellent conver- 
gence of the energy differences, stress tensors and structural 
parameters. Dipole corrections to the electrostatic potentials 
are also included to partly cancel the interactions between the 
slab and its images (221 ESI- The bottom surface is passivated 
by fractional hydrogens I3T1 . The passivated surface is semi- 
conducting, with no surface states in the fundamental band 
gap of GaN. For trial calculations, the slab contains 5 layers 
of GaN, with the top 2 layers being relaxed. To obtain accu- 
rate surface energies, we use a slab containing 10 layers with 
the top 5 layers relaxed. Our tests suggest that the surface en- 
ergy is converged to better than 0.05 eV per surface cell. Our 
surface calculations are restricted to a 2 x 2 or smaller surface 
cell. 



Fig. [4^ gives the convex hull diagram from our EA search, 
and the five corresponding stable structures. Compared to the 
cleaved surface (Fig. [4]:), structure SI has two Ga adlayers. 
Structure S2 has one Ga adlayer. Structure S3 has the top 
N and half of the second N removed. Structure S4 has only 
the top N removed. Structure S5 has an additional N at the 
bridging position of the two top N atoms. The first 4 struc- 
tures have been found by Akiyama et al [32). Structure S5, 
was not reported before and is not intuitively obvious. This 
demonstrates the power of the automated searching by the 
evolutionary algorithm. An analogous Se-trimer feature has 
been predicted to be stable on ZnSe (100) reconstructions at 
Se-rich conditions [ 33 ]. Compared to the cleaved surface, the 
stable structures can be constructed through either vacancies 
or adatoms. No more complicated reconstructions have been 
found. For the Ga-rich conditions, there are many intermedi- 
ate structures between structure SI and S2, lying very close 
to the convex hull (see Fig. [4]). This can be explained by the 
weak directionality of the Ga-Ga bonds. These structures may 
appear under finite temperatures. For the N-rich conditions, 
however, most structures are much higher above the convex 
hull, due to the strong covalent bonding nature of N-N inter- 
actions. 

The features possessed by the stable configurations are: i) 
Ga layers; ii) N vacancies; iii) N trimers. Another common 
feature is N2-dimers l33l (these donot lie on the convex hull, 
but indeed frequently appear as energetically competitive in 
our EA search). 

The electron counting (EC) rule 134U361 emerges from the 
study of semiconductor surfaces. This rule, although not ex- 
act, is very useful in comparing surface energies, estimating 
surface charges and guessing surface reconstructions. It states 
that the surface structures with filled dangling bonds on the 
electronegative elements and empty dangling bonds on the 
electropositive elements are the lowest in energy. The sur- 
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FIG. 4. (a) Energies of GaN(lOll) surface reconstructions explored by the evolutionary algorithm; (b) the corresponding zero temperature 
phase stability diagram up to 2 x 2 supercell reconstructions. Side view and top view of (c) the cleaved GaN(lOll) surface and its lower- 
energy reconstructions found by the EA search; (d) SI configuration: (lxl) Ga-bilayer; (e) S2 configuration: (lxl) Ga-monolayer; (f) S3 
configuration: (2x 1)-1.5 N vacancy; (g) S4 configuration: (2x 1)-1 N vacancy; (h) S5 configuration: (2x1) N-trimer. Large spheres, gallium; 
small spheres, nitrogen. 



face then is likely to be semiconducting. Surfaces not satis- 
fying this rule will be metallic and higher in energy. None 
of the stable reconstructions found here satisfy the electron 
counting rule. Large reconstructions may satisfy the EC rule. 
Indeed, Lee and Kim l37l demonstrated that low-energy semi- 
conducting surfaces can be achieved by large reconstructions, 
complying with EC rule by the modulation of N dimer and 
dimer- vacancy. Due to the limitation of computational power, 
we are only able to predict the reconstructions up to 2 x 2 
surface cells. But our results provide a starting point for pre- 
dicting large reconstructions. 

Let us now employ the same method to search for recon- 
structions of GaN(lOll) surfaces with added oxygen. The 
formation energy can be defined as 

formation = E - /x(Ga)[n(Ga) - n(N)] - /x(0)n(0). (6) 

The oxygen chemical potential /i(O) is a new variable. It has 
no lower bound, but would possibly reach the upper bound by 
the constraints based on the possible O2, N^O^ molecular gas 
and Ga2C>3 solid: 

Mo) < \e{o 2 ), 

2/i(Ga) + 3/i(0) < £(Ga 2 3 ), (7) 
ayi(Ga) + yfi(0) < E(N x O y ). 



Therefore, the convex hull has two dimensions in n(O) and 
[n(Ga) — n(N)]. Fig. [5^ shows the two-dimensional phase 
diagram in space of /i(O) and [/i(N) — /i(Ga)] transformed 
from the convex hull found in our EA search. The two new 
major reconstructions are structure S6 and S7. Compared to 
the cleaved surface, structure S6 has half of the top N re- 
moved, half of the top N and all of the second N replaced 
by O. Structure S7 has the top two N replaced by O. Simi- 
lar reconstructions for the (1011) surface have been reported 
(38). These reconstructions are favored because of the strong 
Ga-0 bond and the electron counting rule. For example, the 
simple cleaved surface has 3 nitrogen dangling bonds per sur- 
face cell. Because each N dangling bond needs 3/4 electron 
to become filled, the surface needs 9/4 electrons per surface 
unit cell to satisfy the electron counting rule. For structure S7, 
substituting 2 N with 2 O introduces 2 extra electrons. There- 
fore, it only needs 1/4 electron to satisfy the electron counting 
rule. 



CONCLUSIONS 

In summary, we present a specially designed method to au- 
tomatically explore low energy surface reconstructions with 
variable surface atoms and reconstruction cells. The power of 
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FIG. 5. (a) Phase diagram of the GaN(lOll) surface with oxygen 
present up to 2 x 2 supercell reconstructions. Important reconstruc- 
tions are labeled from SI to S7. Structures SI to S5 are shown in 
Fig. [4] Side view and top view of structures S6 and S7 are shown 
in (b) and (c). Large balls, gallium; small gray balls, nitrogen; small 
red balls, oxygen. 



this approach is illustrated by the identification of N3 trimeric 
reconstruction of GaN(lOll) surfaces. This method can be 
generally applied to other binary and even multi-component 
systems. 
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